#!/bin/bash

# generates plotting script for spectral density output (Energy components
#
# we assume the different modes of phi and Fphi get dumped to spectra_{nki}_[real/imag].out 
#
# mixed -layer variance goes into ml_var.out
#------------------------------------------------------_#
# globals:
#------------------------------------------------------ #

opt_all_steps=
MLVARdata=ml_var.out
#------------------------------------------------------ #
# read options
while [ "$1" != "" ]; do
  if [ "$1" == "all" ]; then
   opt_all_steps="true"
  fi
  echo $1
  shift 
done;
#------------------------------------------------------ #
gnuscript="plotspec.plt"
gnubin=$(dirname $(which gnuplot))/gnuplot
if [ "$gnubin" == "" ]; then
  echo Error: cannot find gnuplot. > /dev/stderr; exit 1;
fi;

# we also need the transpose, getcolmax, getcolmin awk scripts
getmax=$(which getcolmax)
getmin=$(which getcolmin)
transp=$(which transpose)
qgcmbin=./main_opt
if [ "$getmax" == "" ]; then
  echo Error: getcolmax script must be located in PATH. > /dev/stderr; exit 1;
fi
if [ "$getmin" == "" ]; then
  echo Error: getcolmin script must be located in PATH. > /dev/stderr; exit 1
fi
if [ "$transp" == "" ]; then
  echo Error: transpose script must be located in PATH. > /dev/stderr; exit 1;
fi
if [ ! -f ${qgcmbin} ]; then
  echo Error: $qgcmbin must be located in currenct directory. > /dev/stderr; exit 1;
fi

KMIN=$(${qgcmbin} -eval "K_DOMAIN#1");
KMAX=$(${qgcmbin} -eval "K_DOMAIN#(size(K_DOMAIN))");

echo '#!'${gnubin} > $gnuscript
echo '
# Automatically generated by "'$0'"
set terminal postscript enhanced color lw 2
' >> $gnuscript
##########################################################
# plot spectral density of mixed-layer temp
##########################################################
if [ -f ${MLVARdata} ]; then
 echo '
   set xlabel "k magnitude"
   set ylabel "mixed-layer variance"
      ' >> $gnuscript
   # find appropriate plot range
   numtsteps=$(head -n 2 ./ml_var.out | tail -n 1 | awk '{ print NF-1 }')
   numlines=5
   dTSTEP=$(echo $numtsteps $numlines | awk '{ printf("%i",($1/$2)+1); }');
   echo 'set output "mixed-layer-var.ps"' >> $gnuscript
   echo 'f="./ml_var.out"' >> $gnuscript
   echo 'set title "Mixed-layer isotropic variance spectral decomposition"' >> $gnuscript
   echo -n 'plot ' >> $gnuscript
   tstep=2
   sep=""
   while [ $tstep -le $numtsteps ]; do
    echo -n $sep'f using 1:'$tstep' w lines title "timestep='$(($tstep-1))'"' >> $gnuscript
    tstep=$(($tstep+$dTSTEP));
    sep=","
   done;
   echo -n $sep'f using 1:'$numtsteps' w lines title "timestep='$numtsteps'"' >> $gnuscript
   # do final tstep
   echo >> $gnuscript
fi
##########################################################
# plot spectral density of phi
##########################################################
echo set xlabel \"k magnitude\" >> $gnuscript
echo >> $gnuscript
# plot spectral data
for data in spectra_*_real.out spectra_*_imag.out; do
 if [ -f ${data} ]; then
  noextname=${data%.out};
  nkmode=${noextname%_real}
  nkmode=${nkmode%_imag}
  nkmode=${nkmode#spectra_}
  complextype=${data%.out};
  complextype=${complextype#spectra_${nkmode}_}
  echo set ylabel \"2*PI*k*\sigma\(k\)*Phi^{$nkmode}\" >> $gnuscript
  fname=transposed_${data}
  transpose $data | awk '{ if(NR>1) print $0 }' > $fname
  # do one for each time
  numtimes=`cat $fname | awk '{ x=NF-1; print x; exit; }'`
  echo >> $gnuscript
  ctime=0;
  while [[ $ctime -lt $numtimes ]]; do
   if [[ $ctime -eq $(($numtimes-1)) ]]; then # also plot the overlayed spectrum
    col2=$(head -n 1 $fname | awk '{ print NF }');
    col1=2
    maxvalue1=$(getcolmax $fname $col1 | awk '{ abs=$0; if(abs<0){abs=abs*-1;}; x=$0+0.2*abs; print x; }');
    minvalue1=$(getcolmin $fname $col1 | awk '{ abs=$0; if(abs<0){abs=abs*-1;}; x=$0-0.2*abs; print x; }');
    maxvalue2=$(getcolmax $fname $col2 | awk '{ abs=$0; if(abs<0){abs=abs*-1;}; x=$0+0.2*abs; print x; }');
    minvalue2=$(getcolmin $fname $col2 | awk '{ abs=$0; if(abs<0){abs=abs*-1;}; x=$0-0.2*abs; print x; }');
    maxvalue=$(echo $maxvalue1 $maxvalue2 | awk '{ if($1 > $2){ print $1; }else{ print $2; }}'); 
    minvalue=$(echo $minvalue1 $minvalue2 | awk '{ if($1 < $2){ print $1; }else{ print $2; }}');
    if [[ "$maxvalue" != "$minvalue" ]]; then
      echo set output \"Initial_and_Final_spectrum_${complextype}_nk${nkmode}.ps\" >> $gnuscript
      echo set title \"Spectral Energy density \($complextype\) nki=${nkmode} >> $gnuscript
      #echo set yrange \[ $minvalue : $maxvalue \] >> $gnuscript
      echo "set yrange [ 1e-5 : * ]" >> $gnuscript
      echo plot \"$fname\" using 1:$col2 title \"final tstep\" w boxes, \"$fname\" using 1:$col1 title \"initial tstep\" w boxes >> $gnuscript
      echo >> $gnuscript
    fi
   fi
   if [ "$opt_all_steps" == "true" ]; then
    col=$(($ctime+2));
    maxvalue=$(getcolmax $fname $col | awk '{ abs=$0; if(abs<0){abs=abs*-1;}; x=$0+0.2*abs; print x; }');
    minvalue=$(getcolmin $fname $col | awk '{ abs=$0; if(abs<0){abs=abs*-1;}; x=$0-0.2*abs; print x; }');
    if [ "$minvalue" != "$maxvalue" ]; then
      echo set output \"${noextname}_tstep_${ctime}.ps\" >> $gnuscript
      echo set title \"Spectral Energy density \($complextype\) nki=${nkmode}\, tstep=$ctime\" >>$gnuscript
      # the spectra here are always positive values
      #echo set yrange \[ $minvalue : $maxvalue \] >> $gnuscript
      echo plot \"$fname\" using 1:$col title \"\" w lines >> $gnuscript
    else
     echo Not creating plot $fname because data is all zero.
    fi;
    echo >> $gnuscript
   fi;
   ctime=$(($ctime+1))
  done;
 fi
done;
#-----------------------------------------------------------------------
# now create the tendency spectra
echo set key top left >> $gnuscript
echo set xlabel \"k magnitude\" >> $gnuscript
echo >> $gnuscript
for data in Fspectra_*_real.out Fspectra_*_imag.out; do
 if [ -f ${data} ]; then
  noextname=${data%.out};
  nkmode=${noextname%_real}
  nkmode=${nkmode%_imag}
  nkmode=${nkmode#Fspectra_}
  complextype=${data%.out};
  complextype=${complextype#Fspectra_${nkmode}_}
  echo set ylabel \"2*PI*k*\sigma\(k\)*Phi^{$nkmode}\" >> $gnuscript
  fname=transposed_${data}
  transpose $data | awk '{ if(NR>1) print $0 }' > $fname
  # do one for each time
  numtimes=`cat $fname | awk '{ x=NF-1; print x; exit; }'`
  echo >> $gnuscript
  ctime=0;
  while [[ $ctime -lt $numtimes ]]; do
    if [ "$opt_all_steps" == "true" ]; then
     col=$(($ctime+2));
     maxvalue=$(getcolmax $fname $col | awk '{ abs=$0; if(abs<0){abs=abs*-1;}; x=$0+0.2*abs; print x; }');
     minvalue=$(getcolmin $fname $col | awk '{ abs=$0; if(abs<0){abs=abs*-1;}; x=$0-0.2*abs; print x; }');
     if [ "$minvalue" != "$maxvalue" ]; then
      echo set output \"${noextname}_tstep_${ctime}.ps\" >> $gnuscript
      echo set title \"Spectral Energy density tendency \($complextype\) nki=${nkmode}\, tstep=$ctime\" >>$gnuscript
      #echo set yrange \[ $minvalue : $maxvalue \] >> $gnuscript
      echo plot \"$fname\" using 1:$col title \"\" w lines >> $gnuscript
      echo >> $gnuscript
     else
      echo Not creating plot $fname because data is all zero.
     fi
    fi
    ctime=$(($ctime+1))
  done;
 fi
done;
chmod u+x $gnuscript
./$gnuscript
exit 0;
